import geopandas as gpd
import trimesh
import numpy as np
import os
from shapely.geometry import Polygon as ShapelyPolygon
from shapely.ops import unary_union
from shapely.validation import make_valid
from shapely.errors import TopologicalError

def geojson_to_3d(geojson_path, output_path, scale=100000, export_format='glb'):
    """
    将GeoJSON中的MultiPolygon转换为3D模型，适配新数据格式（用“名称_nam-面积_are”作为小地块命名依据）
    """
    # 检查输入文件是否存在
    if not os.path.exists(geojson_path):
        print(f"错误：输入文件不存在 - {geojson_path}")
        return

    # 创建输出目录
    output_dir = os.path.dirname(output_path)
    if not os.path.exists(output_dir):
        os.makedirs(output_dir)
        print(f"已创建输出目录：{output_dir}")

    # 验证导出格式
    supported_formats = ['glb', 'gltf', 'obj', 'stl', 'ply']
    if export_format not in supported_formats:
        print(f"错误：不支持的导出格式 {export_format}，使用默认格式 glb")
        export_format = 'glb'

    # 确保输出路径包含正确的扩展名
    if not output_path.lower().endswith(f'.{export_format}'):
        output_path = f"{output_path}.{export_format}"
        print(f"调整输出路径为：{output_path}")

    # 检查FCL是否可用
    try:
        import fcl
        fcl_available = True
        print("检测到python-fcl库，将使用精确碰撞检测")
    except ImportError:
        fcl_available = False
        print("未检测到python-fcl库，将使用边界框碰撞检测")

    # ---------------------- 1. 读取并解析GeoJSON（适配新数据字段） ----------------------
    try:
        gdf = gpd.read_file(geojson_path)
        print(f"成功读取GeoJSON：共包含 {len(gdf)} 个小地块")
        # 打印前2个要素的字段信息，确认数据格式
        print("前2个要素的核心字段示例：")
        for idx in range(min(2, len(gdf))):
            row = gdf.iloc[idx]
            print(f"  要素{idx}：名称_nam={row['名称_nam']}, 面积_are={row['面积_are']:.4f}")
    except Exception as e:
        print(f"读取GeoJSON失败：{str(e)}")
        return

    # 检查原始数据重叠
    try:
        union = unary_union(gdf.geometry)
        if union.geom_type in ["MultiPolygon", "GeometryCollection"]:
            print(f"检测到 {len(union.geoms)} 个独立区域，可能存在原始数据重叠")
    except TopologicalError:
        print("警告：原始数据可能存在拓扑错误")

    # 收集所有坐标用于计算偏移量（过滤无效多边形）
    all_coords = []
    valid_geometries = []
    # 必需字段列表（适配新数据：名称_nam、面积_are，用于命名和标识）
    required_fields = ["名称_nam", "面积_are"]

    for idx, row in gdf.iterrows():
        # 跳过空几何
        if row.geometry.is_empty:
            print(f"警告：第 {idx} 个地块为空白几何，已跳过")
            continue

        # 处理MultiPolygon/Polygon类型（新数据geometry可能为MultiPolygon）
        if row.geometry.geom_type == "MultiPolygon":
            polygons = list(row.geometry.geoms)
        else:
            polygons = [row.geometry]

        if not polygons:
            print(f"警告：第 {idx} 个地块没有多边形数据，已跳过")
            continue

        # 检查是否包含所有必需字段（新数据核心字段：名称_nam、面积_are）
        missing_fields = [f for f in required_fields if f not in row or row[f] is None]
        if missing_fields:
            print(f"警告：第 {idx} 个地块缺少必需字段 {missing_fields}，已跳过")
            continue

        # 收集有效坐标（过滤顶点不足的多边形）
        has_valid_coords = False
        for poly in polygons:
            if poly is None or poly.exterior is None:
                continue
            coords = list(poly.exterior.coords)
            if len(coords) >= 4:  # 闭合多边形至少4个点（含起点）
                all_coords.extend(coords)
                has_valid_coords = True

        if has_valid_coords:
            valid_geometries.append(idx)
        else:
            # 生成地块临时名称用于日志（新数据格式：名称_nam-面积_are）
            temp_name = f"{row['名称_nam']}-{row['面积_are']:.4f}"
            print(f"警告：第 {idx} 个地块（{temp_name}）无有效顶点，已跳过")

    if not all_coords:
        print("错误：未找到任何有效的坐标数据，无法创建模型")
        return

    # 计算坐标偏移量（双精度保留经纬度精度）
    all_coords = np.array(all_coords, dtype=np.float64)
    min_lon, min_lat = all_coords.min(axis=0)[:2]
    offset = np.array([min_lon, min_lat], dtype=np.float64)
    print(f"坐标偏移原点：经度={min_lon:.8f}, 纬度={min_lat:.8f}")

    # ---------------------- 2. 逐个构建小地块网格（用新数据字段命名） ----------------------
    scene = trimesh.Scene()
    geometries = []
    existing_meshes = []

    # 检查三角化引擎（确保网格生成稳定性）
    try:
        import mapbox_earcut
        print("将使用 mapbox_earcut 作为三角化引擎")
        os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "earcut"
    except ImportError:
        try:
            import triangle
            print("将使用 triangle 作为三角化引擎")
            os.environ["TRIMESH_TRIANGULATE_ENGINE"] = "triangle"
        except ImportError:
            print("警告：未找到三角化引擎，尝试使用trimesh默认引擎")
            os.environ.pop("TRIMESH_TRIANGULATE_ENGINE", None)

    # 核心配置参数（保持稳定性，适配新数据坐标范围）
    EXTRUSION_HEIGHT = 1e-6  # 极小高度避免平面渲染错误
    SIMPLIFY_TOLERANCE = scale * 0.00001  # 避免过度简化导致形状失真
    MIN_VERTICES_AFTER_SIMPLIFY = 3  # 简化后至少保留3个顶点（构成三角形）
    MIN_DISTANCE_THRESHOLD = 1e-4  # 碰撞检测距离阈值（避免网格重叠）

    for idx in valid_geometries:
        row = gdf.iloc[idx]
        # 提取新数据核心字段生成唯一名称（处理数据类型，确保格式统一）
        地块名称 = str(row["名称_nam"]).strip() if isinstance(row["名称_nam"], (str, int, float)) else "未知名称"
        地块面积 = round(float(row["面积_are"]), 4) if isinstance(row["面积_are"], (str, int, float)) else 0.0
        # 最终命名格式：名称_面积（面积保留4位小数，避免重复）
        网格名称 = f"地块_{地块名称}_面积_{地块面积}"

        # 处理多边形（MultiPolygon/Polygon）
        if row.geometry.geom_type == "MultiPolygon":
            polygons = list(row.geometry.geoms)
        else:
            polygons = [row.geometry]

        if not polygons:
            continue

        # 处理每个子多边形（避免单个地块内多子多边形重名）
        for poly_idx, poly in enumerate(polygons):
            # 严格检查多边形有效性
            if poly is None or poly.exterior is None:
                print(f"警告：{网格名称} 的子多边形 {poly_idx} 无效，已跳过")
                continue

            coords = list(poly.exterior.coords)
            if len(coords) < 4:
                print(f"警告：{网格名称} 的子多边形 {poly_idx} 顶点不足（{len(coords)}个），已跳过")
                continue

            # 修复无效多边形（保留原始形状，避免拓扑错误）
            original_poly = poly
            try:
                if not poly.is_valid:
                    print(f"警告：{网格名称} 的子多边形 {poly_idx} 无效，尝试修复...")
                    poly = make_valid(poly)
                    # 修复后若为MultiPolygon，取面积最大的子多边形（确保主形状保留）
                    if poly.geom_type == "MultiPolygon":
                        largest_poly = max(poly.geoms, key=lambda p: p.area if p.exterior else 0)
                        poly = largest_poly
                    # 重新检查坐标有效性
                    coords = list(poly.exterior.coords) if poly.exterior else []
                    if len(coords) < 4:
                        print(f"警告：{网格名称} 的子多边形 {poly_idx} 修复后顶点不足，已跳过")
                        continue
            except Exception as e:
                print(f"警告：{网格名称} 的子多边形 {poly_idx} 修复失败：{str(e)}，使用原始多边形")
                poly = original_poly
                coords = list(poly.exterior.coords)
                if len(coords) < 4:
                    print(f"警告：{网格名称} 的子多边形 {poly_idx} 顶点不足，已跳过")
                    continue

            # 坐标转换（经纬度→平面坐标，保留双精度避免精度丢失）
            planar_coords = (np.array(coords, dtype=np.float64)[:, :2] - offset) * scale

            # 创建路径并简化（仅顶点数>500时简化，避免小多边形形状失真）
            path = trimesh.load_path(planar_coords)
            original_vertex_count = len(planar_coords)
            if original_vertex_count > 500:
                simplified_path = path.simplify(tolerance=SIMPLIFY_TOLERANCE)
                simplified_count = len(simplified_path.vertices)
                if simplified_count >= MIN_VERTICES_AFTER_SIMPLIFY:
                    path = simplified_path
                    print(f"{网格名称} 的子多边形 {poly_idx} 简化：{original_vertex_count}→{simplified_count} 顶点")
                else:
                    print(f"{网格名称} 的子多边形 {poly_idx} 简化后顶点不足，保留原始顶点")

            # 检查路径有效性（确保可生成网格）
            if len(path.vertices) < MIN_VERTICES_AFTER_SIMPLIFY:
                print(f"警告：{网格名称} 的子多边形 {poly_idx} 顶点不足，已跳过")
                continue

            # 创建网格（挤压生成3D形状，失败时手动兜底）
            try:
                extrusion = path.extrude(height=EXTRUSION_HEIGHT)
            except Exception as e:
                print(f"警告：{网格名称} 的子多边形 {poly_idx} 挤压失败：{str(e)}，手动创建网格")
                # 手动生成三角面（扇形分割，确保平面网格有效）
                vertices = np.column_stack([planar_coords, np.zeros(len(planar_coords))])
                faces = [[0, i, i+1] for i in range(1, len(planar_coords)-1)]
                if len(faces) == 0:
                    print(f"警告：{网格名称} 的子多边形 {poly_idx} 无法创建网格，已跳过")
                    continue
                extrusion = trimesh.Trimesh(vertices=vertices, faces=np.array(faces))

            # 处理网格列表（取第一个有效网格，避免空列表问题）
            if isinstance(extrusion, list):
                if len(extrusion) == 0:
                    print(f"警告：{网格名称} 的子多边形 {poly_idx} 无网格，已跳过")
                    continue
                mesh = extrusion[0]
            else:
                mesh = extrusion

            # 碰撞检测（避免网格重叠，微调位置）
            collision_found = False
            if existing_meshes:
                current_centroid = mesh.centroid
                for existing_mesh in existing_meshes:
                    current_bbox = mesh.bounds
                    existing_bbox = existing_mesh.bounds
                    # 边界框预检查（快速排除无重叠可能的网格）
                    if (current_bbox[0][0] < existing_bbox[1][0] and
                            current_bbox[1][0] > existing_bbox[0][0] and
                            current_bbox[0][1] < existing_bbox[1][1] and
                            current_bbox[1][1] > existing_bbox[0][1]):
                        # 质心距离检查（进一步判断重叠风险）
                        existing_centroid = existing_mesh.centroid
                        distance = np.linalg.norm(current_centroid - existing_centroid)
                        if distance < MIN_DISTANCE_THRESHOLD:
                            # 有FCL库则精确检测，无则直接判定重叠
                            if not fcl_available or (fcl_available and check_collision_fcl(mesh, existing_mesh)):
                                print(f"警告：{网格名称} 的子多边形 {poly_idx} 重叠，微调位置")
                                mesh.apply_translation([1e-5, 1e-5, 0])  # 微小偏移避免重叠
                                collision_found = True
                                break

            # 设置网格名称（子多边形序号确保唯一，避免同地块多子多边形重名）
            mesh_full_name = f"{网格名称}_子多边形_{poly_idx}" if len(polygons) > 1 else 网格名称
            if not hasattr(mesh, 'name'):
                mesh = trimesh.Trimesh(
                    vertices=mesh.vertices,
                    faces=mesh.faces,
                    metadata={'name': mesh_full_name}
                )
            else:
                mesh.name = mesh_full_name

            # 添加到场景和全局列表
            geometries.append(mesh)
            existing_meshes.append(mesh)
            print(f"已创建小地块：{mesh_full_name}（顶点数：{len(coords)}）")

    # ---------------------- 3. 导出3D模型 ----------------------
    try:
        if len(geometries) == 0:
            print("错误：没有生成任何有效的网格，无法导出模型")
            return

        scene.add_geometry(geometries)
        scene.export(file_type=export_format, file_obj=output_path)
        print(f"\n3D模型已导出至：{output_path}")
        print(f"共创建 {len(geometries)} 个小地块网格")
    except Exception as e:
        print(f"导出3D模型失败：{str(e)}")
        # 备选格式兜底（glb失败用obj，反之亦然）
        fallback_format = 'obj' if export_format == 'glb' else 'glb'
        fallback_path = output_path.rsplit('.', 1)[0] + f'.{fallback_format}'
        try:
            scene.export(file_type=fallback_format, file_obj=fallback_path)
            print(f"备选格式导出成功：{fallback_path}")
        except Exception as fe:
            print(f"备选格式导出也失败：{str(fe)}")


def check_collision_fcl(mesh1, mesh2):
    """FCL碰撞检测（有库时使用，精确判断网格是否重叠）"""
    try:
        import fcl
        # 构建FCL碰撞对象（绑定网格顶点和面）
        obj1 = fcl.CollisionObject(fcl.BVHModel(), mesh1.transform)
        obj2 = fcl.CollisionObject(fcl.BVHModel(), mesh2.transform)
        fcl.buildBVH(obj1.collision_geometry, mesh1.vertices, mesh1.faces)
        fcl.buildBVH(obj2.collision_geometry, mesh2.vertices, mesh2.faces)
        # 执行碰撞检测
        req = fcl.CollisionRequest()
        res = fcl.CollisionResult()
        fcl.collide(obj1, obj2, req, res)
        return res.is_collision
    except Exception as e:
        print(f"FCL碰撞检测出错：{str(e)}")
        return False


# ---------------------- 执行转换（适配新数据路径） ----------------------
if __name__ == "__main__":
    # 输入输出路径（根据实际情况修改，新数据为“寒葱2000.geojson”）
    input_geojson = r"C:\Users\ma\Desktop\GIS\寒葱2000.geojson"  # 新数据GeoJSON路径
    output_model = r"C:\Users\ma\Desktop\GIS\寒葱2000"  # 输出3D模型路径（无需加后缀）

    # 执行转换（scale=1e6保留经纬度小数精度，适配新数据坐标范围）
    geojson_to_3d(
        geojson_path=input_geojson,
        output_path=output_model,
        scale=1000000,  # 放大坐标以保留经纬度精度（避免小数位丢失）
        export_format='glb'  # 支持glb/gltf/obj等格式，按需修改
    )